function results_fnames=osl_run_source_recon_inverse(oat)

% results_fnames=osl_run_source_recon_beamform(oat)
%
% takes in an OAT, which needs to be setup by calling oat=osl_setup_oat(S), struct 
% and runs beamformer
% 
% This function should normally be called using osl_run_oat(oat);
%
% MW 2011

global OSLDIR;

source_recon=oat.source_recon;

source_recon.multidipole_method='none';
    
%%%%%%%%%%%%%%%%%
% setup mni coords using std space mask
if(isfield(source_recon,'mask_fname'))

    if ~isempty(source_recon.gridstep)
        mask_fname=osl_resample_nii(source_recon.mask_fname, [source_recon.mask_fname '_' num2str(source_recon.gridstep) 'mm'], source_recon.gridstep, 'nearestneighbour', [OSLDIR '/std_masks/MNI152_T1_' num2str(source_recon.gridstep) 'mm_brain']);
    else
        mask_fname = source_recon.mask_fname;
    end

    disp(['Using mask ' mask_fname]);

    % setup mni coords using std space mask
    [ mni_coord xform ] = osl_mnimask2mnicoords(mask_fname);
    
    str=['Using dipoles at MNI coordinates ' ];

    for vox=1:size(mni_coord,1),
        str=[str ', [' num2str(mni_coord(vox,:)) ']'];
    end;

    disp(str);
elseif(isfield(source_recon,'single_mni_coord'))
    disp(['Using MNI coordinate ' num2str(source_recon.single_mni_coord)]);
    mni_coord=source_recon.single_mni_coord;
    
    mask_fname='';
elseif(isfield(source_recon,'mni_coords'))    
    mni_coord=source_recon.mni_coords;
    
    str=['Using dipoles at MNI coordinates ' ];

    for vox=1:size(mni_coord,1),
        str=[str ', [' num2str(mni_coord(vox,:)) ']'];
    end;

    disp(str);
    
    mask_fname='';    
else
    mask_fname=[OSLDIR '/std_masks/MNI152_T1_' num2str(source_recon.gridstep) 'mm_brain'];
    disp(['Using whole brain']);
    
    % setup mni coords usource_recong std space mask
    [ mni_coord xform ] = osl_mnimask2mnicoords(mask_fname);    
end;

if(~strcmp(mask_fname,''))
   mask=read_avw(mask_fname); 
   save_avw(mask,[oat.source_recon.dirname '/source_recon_mask'], 'f', [source_recon.gridstep,source_recon.gridstep,source_recon.gridstep,1]); 
end;

%%%%%%%%%%%%%%%%%
%% Setup M matrix of MNI coordinate combinations (particulary for
% multi-dipole beamformers)
switch source_recon.multidipole_method  
  case 'none'
    M = mat2cell(mni_coord,ones(1,size(mni_coord,1)),3);
  case 'null'
    % Simplest case - nulling beamformer
    if(size(mni_coord,2)~=3)
        error('mni_coords need to be transposed');
    end;
    if(size(source_recon.null{1},2)~=3 & source_recon.null{1}~=[])
        error('null mni_coord needs to be transposed');
    end;
    M = mat2cell(mni_coord,ones(1,size(mni_coord,1)),3);
    nulls = repmat(source_recon.null,size(M));
    M = [M nulls];
    
  case 'pairwise'
    % All pairwise combinations
    M = mat2cell(mni_coord,ones(1,size(mni_coord,1)),3);
    pairs = combnk(1:numel(M),2);
    M = [M(pairs(:,1)) M(pairs(:,2))];
  case 'l_r'
    % Right-Left hemisphere
    m = [mni_coord,mni_coord*diag([-1,1,1])];
    M = mat2cell(m,ones(1,size(mni_coord,1)),[3 3]);
    tol = 10; %mm
    M(abs(m(:,1)) < tol,2) = {[]};

    source_recon.beamformer_fn=@osl_beamformer_fuse;
  case 'predefined'
    if isfield(source_recon,'M')
      M = source_recon.M;
    else
      error('You must define M matrix if method set to ''predefined''')
    end
end

if(size(M,1)<5 && size(M,2)<5)
    cell2mat(M)
else
    sizeM=size(M);
    sizeM
end;

for subi_todo=1:length(source_recon.sessions_to_do),   
        
    subi=source_recon.sessions_to_do(subi_todo);

    source_recon.session_name=['session' num2str(subi)];
        
    disp(['Using ' source_recon.method ' for source recontruction']);  % changed dy DM

    source_recon_sub=source_recon;
    source_recon_sub.do_plots=oat.do_plots;

    if ~isempty(source_recon.D_continuous),
        source_recon_sub.D_continuous=source_recon.D_continuous{subi};
    else
        source_recon_sub.D_continuous=[];
    end;
    
    if ~isempty(source_recon.D_epoched),
        source_recon_sub.D_epoched=source_recon.D_epoched{subi};
    else
        source_recon_sub.D_epoched=[];
    end;
    
    source_recon_sub.mri=source_recon.mri{subi};

    source_recon_sub.mni_coord=cell2mat(M);

    if(isfield(source_recon,'hmm_block'))
        source_recon_sub.hmm_block = source_recon.hmm_block{subi};
    end;
    
    source_recon_results=osl_inverse_batch(source_recon_sub);  % changed by DM
    %source_recon_results=feval(source_recon.beamformer_fn,source_recon_sub); 

    source_recon_results.source_recon=source_recon;
    source_recon_results.mask_fname=mask_fname;
    source_recon_results.mni_coord=mni_coord;        
    source_recon_results.gridstep=source_recon.gridstep;        
    source_recon_results.recon_method=source_recon.method;

    % save

    source_recon_results.session_name=source_recon.session_name;    
    source_recon_results.fname=[source_recon_results.session_name '_recon' ];
    disp(['Saving source-space results: ' source_recon_results.fname]);  % changed by DM
    
    osl_save_oat_results(oat,source_recon_results);

    results_fnames{subi}=source_recon_results.fname;
    
end;
